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O-t- ABSTRACT 

2 , Aims. We calculate the hydrogen and helium ionization in B[e] envelopes and explore their dependence on mass-loss and 

■ effective temperature. We also present simulated observations of the Ha emission line and the CIV AA1550 doublet, 

^ ' and study their behavior. This paper reports our first results in an ongoing study of B[e] supergiants, and provides a 

glimpse on the ionization of the most important elements in self-consistent numerical simulations. 
Methods. Our newly developed 2D stellar atmosphere code, ASTAROTH, was used for the numerical simulations. 
The code self-consistently solves for the continuum radiation, non-LTE level populations, and electron temperature in 
^ ' axi-symmetric stellar envelopes. Observed profiles were calculated by an auxiliary program developed separately from 

O ! ASTAROTH. 

, Results. In all but one of our models, H remained fully ionized — only for M>10^^ Mq yr~^ and Teff <18,000 K did 

OO ' we obtain a neutral H disk, and then only for radii beyond 3R, . Due to ionizations from excited states it is much more 

' difficult to get a H neutral disk than indicated by previous analytical calculations. Near the poles, the ionization is high 

in all models, while helium recombined in the equatorial regions for all but our lowest mass-loss rate (10~® M© yr~^). 
Although the model parameters were not adjusted to provide fits to any particular star, the theoretical profiles show 
some features seen in the profiles of R126. These include the partially resolved double peaked profile of Ha, and the 
weak emission associated with the UV C iv resonance line. 
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H 1. Introduction firmed by interferometr ic observations (jStee et al.l 119951 : 
■ - '• iQuirrenbach et al] |1997[ ). Similarly, a growing body of 
Massive stars are very important constituents of the evidence suggests that the LBV phenomenon includes 
Universe despite their low share of galactic masses. They 2D or 3D processes . For example, tj Car is a binary 
are the primary sources of elements heavier than Li and (iPamineli et alllSl) and the Homunculus nebula around 
power the internal evolution of galaxies by their radiation, ^^e most cited example of a bipolar outflow. The as- 
winds, and explosions. Many of these stars, and their im- pherical envelope of -q Car is not unique among LBVs; 
mediate environment, cannot be fully understood in the |Groh et al.l ^M) found that AG Car is a fast rotator, and 
context of plane-parallel or spherical models, and their mul- ^ence may have a latitude dependent wind. The relevance 
tidimensional nature has to be taken into account if their of the binarity and fast rotation in the LBV phenomenon 
evolution or physical state is to be adequately described. general is not yet clear and is the subject of vig orous 
The origin of their asphericity is either fast rotation, the research (e.g.. iMartin [25061 : [Nielsen et al.l[200l . 
presence of a dynamically important magnetic field, or their g^gj s tars, or rathe r the stars that show B[e] phe- 
mteraction with a companion. nomenon (iLamers et al.lll9il . are characterized by strong 
Three groups of such stars are particularly relevant for broad H Balmer emission lines (as are the classical Be 
this paper. These are the classical Be stars, the Luminous stars) and by the presence of narrow permitted and forbid- 
Blue Variables (LBV), and the B[e] supergiants (sgB[e]), den emission hues from low-ionization species, like Fell, 
all of which share many similar characteristics. The 2D [peH], or [01]. B[e] stars also show strong near/mid-IR ex- 
or 3D nature of these objects was recognized through cess that is attributed to hot cirumstcl lar dust. Intereste d 
both observations and theoretical calculations. For ex- readers should refer t o [zTckgraf et all ([19851 [1981 Il992f ): 
ample, the presence of disks around classical Be stars iMiroshnichenko et all (|2005() and references therein for 
was inferred from line modeling a nd polarimetric stud- further information on the B[e] phenomenon. A particu- 
ies ([Poeckert fc Marlborougg|1978a[|g), and has been con- lar ly interesting subclass of the B[e] stars are the sgB[e]- 
s (ILamers et al.l il998i) . These are single B supergiants 
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|1984[) . Their location on the HR diagram may suggest a 
link to the more massive LBVs or, alternatively, they may 
represent a second e volutionary path between Of and Wol f- 
Rayet (W-R) stars (|Schulte-Ladbecklll998HZickgrajll998^ . 

The envelope of sgB[e] stars is thought to be lat- 
itude dependent, with a normal B supergiant wind 
at the po le and a deii s e an d slo\Y ly moving equato- 
rial flow (IZickgraf et al.l 1985 . Il986l). Polarimet r ic ob - 
servations bv iMagalhaesI (|1992f ): lOudmaiier et all (|l998t) : 
iMelgarejo et al.l ( 200lf ) revealed that sgB[e]-s have strong 
intrinsic polarization consistent with this picture. The most 
widely accepted theory for the bi-modal nature of the enve- 
lope is the rotationally induced bi-stability. Around T^ff— 
20,000 - 25,000 K the characteristics of the line-driven 
wind a bruptly change due to the recombination of Fe IV t o 
Felll (iPauldrach fc Pulsl[T99fll : IL amers fc Pauldrachll991h . 
If the sgB[e] stars are fast rotators, as suggested, the en- 
suing gravity-darkening can place the equatorial and po- 
lar regions on the opposite sides of the bistability, hence 
providing a mechanism to produce the bi-modal envelope. 
However, the bi-stability model is not without its weak- 
nesses. For example, gravity darkening reduces the flux 
at the equator, potentially inhibiting any equatorial en- 
hanced mass-flux. It is also difficult to maintain multi- 
ple scattering in the equator when the photons can eas- 
ily escape in the polar direction, (see discussion in fJ7]). 
Consequently, other sgB[e] models have been proposed, like 
the presence of a Keplerian disk similar to those around 
classical Be stars. Fur ther, disk winds as proposed by 
lOudmai ier et al. (19 98|), the traditional wi nd compressed 
disk model of Bio rkman fc Cassinellil ( 1993f). and even the 
magnetically confined disk mo dels ( ud-Doula fc OwockH 
l2002t lOwocki fc ud-Doulal l2004 ) cannot be excluded. So 
far, theoretical efforts have primarily focused on under- 
standing the hydrodynamic structure in the context of 
the bi-stability model , and trying to inc orporate the bi- 
stability jump into thelCastor et al] (^1975^ formalism (e.g., 
iPelupessv etall 120001 : ICure et al.ll2005D . 

There have only been a few efforts to understand the 
ionization structure of the wind and to perform spectral 
analysis of B[e] stars. Semi-analytical calc ulations of hydro- 
gen an d he lium io n izatio n were done by iKraus fc Lamer j 
(|2003D and iKraud (|2006D . and they found the equatorial 
disk essentially neutra l in the s e spe cies all the way down 
to the stellar surface. I Ported (|2003l ) attempted to repro- 
duce the optical to infra-red continuum of R 126 (prototype 
sgB[e] in the LMC) by using bi-stability and Keplerian vis- 
cous disk models. Both of them reproduced the optical and 
near-IR emission, but they failed to account for the dust 
emissi on by an order of magnitude. Recently, iKraus et al.l 
(|2007t ) proposed that the dust and optical/near-IR contin- 
uum emission can be reconciled if the free-free and bound- 
free emission originates from the polar wind rather than 
the equatorial disk. 

A self-consistent spectral analysis of a stellar envelope 
needs to solve the coupled equations for the radiation field 
and level populations. Obviously, this is a very difficult 
task in a non-LTE 2D or 3D en velope and compromises 
and simplif i cation s are inevitable. iKraus fc Lameri (|2003l ) 
and iKraui (|2006[ ). for example, used a simplified radiative 
transfer, limited their analysis to optically thin or optically 
thick cases, and neglected ionization from excited levels; all 



gained by simple calculations, one needs numerical simula- 
tions for a fully self-consistent spectral analysis. 

To provide a tool for such studies of Be/B[e] stars, bi- 
naries, and LBVs, we have been developing a code for axi- 
symmetric models. So f ar, only basic tests and code verifi- 
catio n were performed (jGeorgiev et al.|[2006t IZsargo et alj 
'2006!) , but we have now reached the point where scientif- 
ically meaningful simulations can be performed. The ion- 
ization structure of hydrogen and helium in axi-symmetric 
sgB[e] envelopes offers such an opportunity. T h e deta ils of 
the prograni are d iscussed in iGeorgiev et al.l (|2006( ) and 
'Zs argo et al.l (j2006l ) so we only briefly describe it in The 
hydrodynamic structures and the atomic models that are 
used in our simulations are discussed in computational 
issues are discussed in !j4l and our results are presented in 
[JSl Simulated observations are shown in fjSl and we draw 
our conclusions in fJTjand fJ51 

2. The Code 

Our C++ code, ASTAROTH, was developed for simula- 
tions of stellar envelopes. However, it is flexible enough to 
be applied to any hot axi-symmetrical object with velocity 
gradients (including extragalactic objects, like AGN-s). 

The non-LTE level populations, the radiation field, and 
the electron temperature are calculated by simultaneously 
solving the equations of statistical equilibrium, radiative 
transfer, and ene rgy conservation. The short-characteristic 
method (see e.g., iMihalas et ai][T978H K unasz fc Aueil[T98i 



iBusche fc Hillied 12000') is used to treat the continuum ra- 
diation transfer while bound-bound transitions are treated, 
for simplicity, by the Sobolev approximation. The simul- 
taneous solution of the equations is found by an a pprox - 
imate lambda iteration (|Rvbicki fc Hummed 119911 Il992t ). 
The code has been tested by solving 2D pure scattering 
problems with grey opacity, as well as reproducing spher i- 
cal symmetric models of CMFGEN (jHillier fc Milledll998D . 
a well-established code in stellar studies. In these tests, the 
new code reproduced the reference results within a few per- 
cent f see. IZsargo et al.ll2006l: IGeorgiev et al.ll2006[ ). 

3. Models 

To simulate the hydrodynamical st ructure of a sgB[e ] atmo - 
sphere we followed the approach of lKraus fc Lamera ()2003l ). 
but relaxed two of their simplifications. We used a /3-law 
(jCastor et al.lll975l ) to describe the radial velocity and al- 
lowed for a varying electron temperature. Our velocity field 
was still simplified; only a latitude dependent radial velocity 
was included and the azimuthal and latitudinal velocities 
were set to zero. To simulate the bi-modal wind, we de- 
scribed the radial velocity and the mass-loss per unit solid 
angle by 



Vr{R, 

and 



-ft* 

1 — - 

R 



(0 = 0)- 10 
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respectively. The parameters R and 9 are the traditional 
polar coordinates, and s controls the thickness of the equa- 
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Table 1. Description of the Models 
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Fig. 1. The density as a function of latitude at R— 100 i?* 
in our models A and B (see Table [J). Note that the density 
scale is logarithmic. The density distribution for models 
C and D is similar in shape, but the values are 10 times 
greater. 
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Fig. 2. A meridional snapshot of the logarithmic hydrogen 
density near the star for models A and B. The pole is to- 
ward the top of the page and the equator is horizontal. The 
envelope is axi-symmetric around the pole and top-bottom 
symmetric around the equator. The units are in stellar ra- 
dius and the thick circle at the center represents the surface 
of the star. A similar snapshot for models C and D would 
look identical except with larger density values. Note, that 
the open contour lines near the stellar surface are artifacts 
of omitting the region R<1.1 from the plot which was 
done for clarity. 



where 10 and 100 were used. In these initial calculations 
we avoided the higher values because of the need to have a 
much finer spatial grid which would have substantially in- 
creased the computational effort. In a future, detailed anal- 
ysis of sgB[c]s, the effects of varying disk thickness will need 
to be explored. 

Using Eq. [5] and assuming top-bottom symmetry, the 
total mass-loss rate can be calculated by 



M = 47r ■ 



dtdn 



(0) ■ sine de 



(3) 



Also, because all non-radial velocities are zero, the gas den- 
sity in the wind, p{r, 6) is given by 



(0) 



(4) 



In our models we used b^ = — 2 and hm— 1 which resulted 
in a 3 order of magnitude density enhancement and a 2 
order of magnitude velocity decrease between the pole and 
the disk (see [JT] for discussion on the validity of Sobolev 
approximation in the disk). The 2D density structure is 
displayed in Fig [1] and Fig O The bi-modal structure of 
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' - Power for the standard velocity law (jCastor et al.lll975f ). 



There i s a significant dif ference between our approach 
and that of Kraus fc Lame rs (2003^. They did not deal with 
the photosphere, and used .Kurucd (|1979D spectra to sim- 
ulate the stellar radiation field. ASTAROTH, similarly to 
CMFGEN, solves for the unknown level populations and 
ionization structure, and calculates spectra for both the 
photospher e and the wind. It u ses the same methods as 
CMFGEN (|Hillier fc Milleilll998t) to create an approximate 
hydrostatic photosph ere. In brief, the win d follows a simple 
beta-like velocity law ijCastor et aLlllQTSl ) , and for this work 
starts at a radius Rptot defined by Vr{Rphot, &)— I km s^^. 
Below this radius the density follows an exponential law 
and the velocity is set by the continuity equation for the 
given mass-loss rate. The value of Rptot is latitude depen- 
dent in the models presented here, but it was always nearly 
i?, (<l-2% difference). 

For the stellar parameters, listed in Table [U we chose 
values that are bro adly representative o f the most lumi- 
nous sg B[e]-s (see, IZickgraf et all Il985l I1986D and some 
less luminou s LBVs (e.g., P Cveni. lPauldrach fc Puls|[l990l : 
iDrewl [l985h . The mass -loss range cover s the u pper end 
of the range used in iKraus fc Lamer j ()2003[ ) and in- 
cludes their baseline models (models A, C, and F in 
Kr aus fc Lame rs 2003). T he suspected ma ss-loss rate of 
R 126 (~4xl0-^ M0 vr- MBiorkmanlll998f) is also within 
our mass-loss range. We vary the effective temperature be- 
tween Teff ~22,500 K (L= 1.5x10^ Lq), that of R 126, to 
Teff 18,000 K (L=6xl0^ L0) which approximately cov- 
ers t he sgB[e] temperatu r e ran ge in the Magellanic Clouds 
(see, IZickgraf et al][l985l . fT986h . 

Finally, the atomic model is presented in Table [21 We 
used essentially H/He atmospheres with a fairly large num- 
ber of levels included. A few levels and ionization states of 
carbon were also included to achieve more realistic heating 
and cooling terms. The fractional (number) abundances of 
He and C was 0.1 and 9.82x10"^*, respectively. The later 
number represents the net abundance of CNO elements. 



4. Computational Issues 

Previous works on st ellar winds and sgB[e] stars (e.g., 
iKraus fc Lamer j|2003[ ) have found that ionization changes 
can occur suddenly, in the form of ionization fronts. One 
of our main concerns was, therefore, whether our spatial 
grid would adequately resolve a potential 2D ionization 



4 



Zsargo et al.: Ionization in the wind of sgB[e] stars 



Table 2. Atomic Model 



Specie 


Number of Levels 


HI 


20 


HII 


1 


He I 


11 


Hell 


20 


He III 


1 


CII 


9 


cm 


10 


CIV 


5 


cv 


1 



semi-irregularly with the limitation that all latitude must 
have the same radial grid. The only problem we encountered 
was the sharp hydrogen recombination front that occurred 
in model D, and the spatial grid had to be adjusted by hand 
to better resolve the front (adaptive gridding capability is 
not yet included in ASTAROTH). Otherwise the standard 
grid was adequate for our models. Simulations on a denser 
grid with 120 depth and 20 latitude points revealed no qual- 
itative changes in the populations and temperatures. 

The convergence was also sometimes slowed by oscil- 
lations in the low radiation field/high density equatorial 
disk. As the statistical equilibrium and energy conserva- 
tion equations are highly nonlinear in the populations, 
radiation field, and electron temperature, solution tech- 
niques can provide oscillating solutions. Often the oscil- 
lations arose at depths, and in populations, where changes 
would have a negligible influence on the emergent spectrum. 
Unfortunately, since we currently use a stopping criterion 
based on the maximum population change, such oscilla- 
tions can greatly increase the total computational effort. To 
dampen the oscillations we linearized the statistical equi- 
librium equations with populations averaged over several 
previous iteration cycles. Switching to linear interpolations 
of quantities, instead of the standard cubic or parabolic 
approximations, also improved the convergence. The lin- 
ear interpolation is an extremely well-behaved and stable 
approximation and it would be the preferred method if it 
provided the necessary accuracy. Unfortunately, this is not 
the case for spatial grids that are practical for ASTAROTH 
simulations; therefore, linear interpolation was limited to 
grid-points immediately around the trouble spot. 

5. Ionization and Temperature Structure 

Figs. El HI and [5] show the electron temperature, hydrogen, 
and helium ionization structures, respectively, in our model 
envelopes. 

5.1. Electron Temperature 

The electron temperature varies between ^100,000 K to 
^^20,000 K in the inner hydrostatic atmosphere (not shown 
in Fig. [3]), and it was a few 1,000 K in the outer envelope. 
The behavior of Tg in the wind is similar, as expected, for 
all models — the equator is cooler than the pole at the 
same radius. The effect of the lower luminosity in models B 
and D is a shift of the overall structure closer to the stellar 
surface (see Fig [3]), while increasing the mass- loss results 
in a thicker cool disk. It is immediately obvious from Fig [3] 



The temperature in our models is too hot for dust for- 
mation. Nowhere in our models the temperature falls bc- 
low ~ 1,500 K which is the upper limit for dust formation 
()Porterj|2003f) . We cannot conclude, however, that dust is 
formed beyond R > 100 — 200i?* (the outer boundary of 
our models) because of the relative simplicity of our models 
and the neglect of adiabatic cooling in our simulations. We 
will address this question in follow-up studies with models 
constrained by the observations of individual stars. 

5.2. Hydrogen Ionization 

Fig. [4] shows the ionized to neutral H ratio in logarithmic 
scale where neutral H means the total population of all 
H I levels. It was very difficult to display the wide range of 
ionization levels occurring in our models. Since there were 
orders of magnitude differences even within a single model, 
we opted for using a separate set of contour levels for each 
model. 

The hydrogen ionization shows latitudinal variations 
similar to those of the temperature. The equatorial region 
is more neutral than the pole in all models, and the lower 
the effective temperature or higher the mass-loss the more 
neutral the inner envelope. In only one model, model D, 
does a neutral H disk form, and then only beyond ~ 3-R^,. 
This contradicts the result of iKraus" fc LamersI (|2003[ ) who 
found predominantly neutral hydrogen disks, even nearly 
at the surface, for very similar models. In our models A, B, 
and C the neutral hydrogen is negligible compared to H II 
(inside lOi?,), although there is a tendency for the neutral 
hydrogen fraction to rise at larger radii. The high level of 
ionization prevails even if the effective temperature is low- 
ered (model B) or the mass- loss rate is increased (model C). 
Only the combined effect of the two (model D) produced 
a neutral hydrogen disk. We will furthe r discuss the dif- 
ferenc es between our results and those of lKraus fc Lameri 
(I2003D in 

5.3. Helium Ionization 

Fig. [5] shows the logarithm of the ionized to neutral He ra- 
tio. Ionized He means the total population of all Hell and 
He III levels. The ionization level of helium is markedly 
lower than that of hydrogen, a result to be expected given 
the low effective temperatures which were adopted. As op- 
posed to H, He forms a neutral disk in all but one of our 
models — in model A, only an enhancement of neutral he- 
lium occurs in the equatorial region. The neutral disks show 
characteristics similar to those of the temperature distribu- 
tion or hydrogen ionization; e.g, lowering Tg/ / results in an 
inward shift of the neutral regions as seen in the last two 
panels of Fig. [5] Higher mass-loss causes the neutral disk 
to become thicker. It also appears that the thicker the disk, 
the sharper its ionization boundaries. The neutral He disk 
reaches down deeper than the hydrogen disk in model D, 
but it is still truncated at ~2 i?*. 

6. Observed Profiles 

Observed profiles for a converged ASTAROTH model are 
calculated independently by an auxiliary routine (described 
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Fig. 3. The temperature structure near the star in models 
A, B, C, and D (top to bottom). The layout of the plots 
is the same as that of Fig. [2l The contour levels are from 
6,000 K to 16,000 K in 2,000 K increments and shown in 
10,000 K units. 



the results of the auxiliary routine. Below we discuss the 
properties of the observed Ha profile, which arises from the 
disk, and the C IV AA1550 doublet profile, which primarily 
arises in the polar wind. We also discuss similarities and 
differences between the computed lines, and those seen for 
our benchmark B[e] supergiant R 126, although we stress 
that this work is not a spectral analysis of this star. For 
the calculations presented in the following sections we gen- 
erally adopted lOkms^^ for the Doppler parameter of the 
intrinsic line absorption/emission profile. Because of the 
low velocities in the disk, the adopted Doppler parameter 
can have a substantial influence on the line profile shape, 
and model profiles calculated with Sobolev approximation 
are significantly different. Furthermore, as the azimuthal 
velocities are zero in our models, the profiles are not rota- 
tionally broadened. 

6.1. H Line Profiles 

The shape of the Ha lines in Fig. [S] are very sensitive to the 
observer's viewing inclination (i). Depending on the model 
either a single or double emission peak is seen; sometimes 
the emission peaks are well resolved, while at other times 
they blend into a continuous emission feature. For i = 0, 
the half-widths at the line base (HWLB) are 30-50 km s'-^ 
which is slightly larger than T^o in the equatorial plane 
(20 km s~^). Very weak wings, arising from emission in 
the wind outside of the equatorial latitudes, is also present. 
At higher inclinations the profiles are broader, and gen- 
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Fig. 4. The log ^ ) contours for models A, B, C, and 

D (top to bottom) . The layout of the figures is the same as 
that of Fig. [2] The thin solid lines represent ratios above 
(ionized H) and the thick dotted lines show neutral regions. 
Note, that the lower the ratio the more neutral the hydro- 
gen is! The spacing between levels is different from figure 
to figure and between neutral and ionized regions. 
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tic feature often seen in r eal sgB[e] observations (see e.g., 
iZickgraf et allllQSSL fl986h . 

The double-peak profiles arise from optical depth effects 
in our simulations, and are not the result of rotation. The 
keys to understanding the double-peak profile are (1) that 
Ha is optically thick at the line center in all four mod- 
els, and hence photons generally escape in the line wings 
where r ~ 1, and (2) that the Doppler velocity is compara- 
ble to the disk outflow velocity. Let's consider a spherical 
model with an extended atmosphere but no velocity field. In 
general, the source function will decrease with radius, and 
thus the maximum of the line source function, along a line 
of sight parallel to z-axis, will be in the xy plane. Optical 
depth unity will occur in the xy plane at two wavelengths, 
symmetrically located about line center. This will produce 
a distinct double-peaked profile with, in the absence of a 
velocity field, both peaks being of equal strength. However, 
when a velocity field is present, the red peak will be stronger 
than the blue. Photons emitted on the red side only see ma- 
terial moving away, and never come into resonance. On the 
other hand, photons emitted in the blue side must traverse 
material that is resonant with it, thus reducing the line in- 
tensity. The separation of the two peaks depends on the 
line width, and the overall profile and equivalent width is 
very sensitive to the adopted Doppler and wind velocities. 
When the wind velocity is much larger than the Doppler 
velocity we see a classic wind emission line with an optical- 
dependent redshift of 1 to 3 Doppler widths. As apparent 
from the above discussion, the blue/red peak asymmetry is 
not a specific feature of 2D models and also works for spher- 
ical envelopes, provided that the line is optically thick and 
the Doppler and flow velocities are comparable. 

The shape of the pole-on viewed Ha in our model A is 
very similar to that s een in the spectra of R 126 (see Fig. 3 of 
IZickgraf et al.|[r985D . However, the simulated profile lacks 
the broad wings of the observations and is a factor of ~10 
weaker than it should be. This indicates that the density is 
much lower in Model A than it is in the envelope of R 126, 
and highlights the importance of frequency redistribution 
in electron scattering (the profiles in Fig. [H] were calculated 
ignoring electron-scattered line photons) . The higher mass- 
loss rate in models C and D increases the equivalent width 
(it still falls short by a factor of ~4, except at i= 90°), 
but the profiles now are very different from that of R 126. 
The very strong absorption around and blueward of the 
line center is not observed. This central absorption is very 
prominent for models B, C, D, and may indicate that hy- 
drogen is more neutral in the line forming region of these 
models than it is around R 126. 

Table [3] lists the equivalent widths (EW), the contin- 
uum levels across Ha, the total intensities arising in the 
Ha profiles (line fluxes), and the Ha to H/3 line flux ratios. 
The equivalent width varies with viewing angle, although 
the variation is weak for low to moderate inclinations. The 
continuum flux also varies with inclination angle. In mod- 
els D and C the variation in continuum flux is large (the 
continuum is '^lOx smaller for the i — 90° model than it 
is for i = 45°), and this causes the EW to be higher for the 
i = 90° models, despite the lower line flux. 

In modeling of O stars, LBVs, and W-R stars simple 
scaling laws can provide valuable insights into wind spec- 
tra, and greatly facilitate modeling through a reduction in 
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Fig. 6. Simulated observations of Ha for models A, B, C, 
and D (top to bottom). There are three viewing angles, i= 
0° (pole-on, sohd fine), 45° (dashed line), and 90° (edge-on, 
dash-dotted line). Note that all profiles are rectified. 



same effective temperature, and the same "optical depth 
invariant" parameter, Q = -, — ^-tt-s-, showed similar spec- 
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Table 3. Ha Equivalent Widths and Line Fluxes 



B 



C 



D 



i ' 




EW (A) 


Fi = EWxFe 


1 1 (H a ) 


0° 

4o 

90° 


1060 

8/0 
384 


61 
oo 
21 


64681 
4or36 
8198 


3.8 

Q O 

o.o 
5.3 


0° 
45° 
90° 


754 
615 
267 


54 
49 

30 


40754 
30320 
8013 


4.3 
4.4 
7.1 


0° 
45° 
90° 


3177 
2263 
107 


195 
204 
873 


618689 
461810 
93375 


5.7 
5.7 
6.4 


0° 
45° 
90° 


1725 
1195 
134 


116 
134 
392 


200859 
160345 
52511 


6.4 
6.6 
7.1 



Note that the values (except the ratios) in this table are 

rounded up to the closest integers! 

^ : Viewing angle, i= 0° (pole-on), 45°, 90° (edge-on). 

: Continuum intensity across Ha (in Jansky). Constant to 1- 
2 %. 



low because the stellar radiation is seen through the dense 
equatorial disk and the projected area of the disk is small. 
The very large EWs for i— 90° viewing angle (see Table [3]) 
are principally the result of the disproportionately low con- 
tinuum; for high inclination angles the continuum weakens 
more with inclination angle than does the Ha emission. 
Since the continuum is low at every wavelength, these stars 
may appear unusually faint, and thus their luminosities and 
mass are underestimated. 

A further important feature of the Ha emission is its 
Teff dependence which mainly affects the absorption at 
and blueward of the line center. The lower the effective 
temperature is the stronger the absorption (see, for exam- 
ple, the difference between models B and A in Fig. [S]). This 
is also apparent on the lower equivalent widths for models 
B or D relative to their respective high Teff counterparts. 
It is not surprising that the hydrogen in these models is also 
less ionized and excited (relative to the models with higher 
Teff) which is clearly shown by the Ha to H/3 flux ratios 
listed in Table [31 The larger population in the low-lying 
states can also increase the self-absorption in the Balmer 
lines; hence it can explain the stronger absorption in models 
B and D. 



Il998( ). With the added complexity of 2D models, the iden- 
tification of useful invariants is not so easy. 

The extra difficulty of the 2D models is apparent from 
Table [21 First, continuum fluxes, line fluxes, and line EWs 
vary with the observer's inclination. Thus knowing the ob- 
server's inclination is crucial if meaningful envelope param- 
eters are to be derived from the spectra. It is also apparent 
that for B[e] stars ionization effects are important, even for 
hydrogen. Models C and D, for example, have identical Qs 
but different EWs and line fluxes. Another complexity, not 
seen in W-R models, is that the EW and line-flux of Ha 
depends on the adopted Doppler velocity. When spectra for 
the very same models but with Vturb= 5 km were cal- 
culated (not shown), the EWs and line fluxes decreased by 
~40%. The Ha to H/3 ratio, on the other hand, was not 
affected by the lower turbulent velocity. 

To fully explore the effects of parameters on the observ- 
ables is beyond the scope of this paper but simple trends 
can already be identified. It is striking, for example, that 
the ratio of the continuum intensities for i— 45° and 0°, 
especially in models C and D, is close to cos(45°)= 0.71 
which is also the ratio of the projected area of a disk to 
its total area when tilted by 45°. One would expect such 
ratios if a substantial portion of the continuum around Ha 
formed in the transition region between the polar wind and 
the dense disk. The visible area of this disk "surface" would 
then control the continuum level. Furthermore, most of the 
Ha emission has to originate from the equatorial region 
since many of its characteristics cannot be explained in the 
context of a nearly spherical polar wind; e. g., its sensi- 
tivity to the observer's inclination angle. Simulations with 
CMFGEN (performed by the authors) also showed that one 
would need extreme stellar parameters to reproduce ob- 
served sgB[e] features in spherical models. Mass- loss rates 
in excess of 10"'' Mq yr~^ and Voo < 200 km s"^ are re- 
quired to match the observed linewidth and EW (~900 A 
■ IzTckgraf et al..,1985.) for R 126. 



6.2. Civ Profile 

The C IV doublet profile, for models A, B, C, & D, and for 
i=0°, 45°, & 90° is shown in Fig. [71 A striking feature of 
the observed profiles in all models is the lack of red-shifted 
emission at pole-on and intermediate viewing angles (i.e., 
1=0° and 45°). This is also one of the notable features in 
the International Ultraviolet Explorer (lUE) observations 
of R 126 (Fig. 5 in lZickgraf et~aLlll985D . 

The behavior of C IV is the exact opposite of what we 
saw in Ha. C IV is mainly in absorption formed in the polar 
wind, and shows little sensitivity to the viewing angle (at 
least until very nearly edge-on view). This is because the 
polar wind occupies a large solid angle around the star and 
it is also nearly spherical. A viewing angle of 0° or 45° 
makes little difference, apart from the position of the blue 
absorption edge. The role of the dense equatorial region in 
the formation of the C IV lines is to block the emission from 
the far side of the envelope. The profiles in Fig. [3 are also 
very sensitive to the effective temperature of the models. 
The absorption is much weaker for the lower Teff models 
(B and D) than for their high temperature counterparts (A 
and C). 

The edge-on viewed spectra in Fig. [71 are again complex. 
There is always a broad and strong emission component, 
the signature of the polar wind. They appear to be very 
strong in the rectified spectra only because they are seen 
against a very weak continuum. In three of the predicted 
profiles, doublet absorption components can also be seen. 
Their strength and shape is highly model dependent, and 
presumably arises from "disk" absorption of both line and 
continuum photons. The C IV profiles potentially offer an 
opportunity to constrain the thickness of the equatorial disk 
if the star is viewed nearly edge-on. 

7. Discussion 
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Fig. 7. Same as Fig.lHbut for CIV AA1550. Note that the 
velocity scale is centered on the doublet line with shorter 
wavelength. 



ex planation. We think that the nebular approximation used 
by iKraus fc Lameri (j2003| ). is the reason for the differing 
results. In the nebular approximation it is assumed that all 
recombination to n> 1 levels of hydrogen will cascade down 
to level n=l and the photo- ionizations from the n>l levels 



Fig. 8. The ratio of the total H photo-ionization (PR) from 
levels with n>l to the total recombination (RR) to levels 
with n>l for the pole (solid line) and equator (dashed line) 
of model A. The thick lines are for ASTAROTH (the actual 
2D model) and the thin curves are the CMFGEN results 
(spherical models created for the given latitude). 



and O stars. iDrewl (|1985[ ) studied P Cy gni, an LBV with a 
mass-loss rate and luminosity similar to those of our mod- 
els, and found that photo-ionization from the n>l states 
were crucial in determining the hydrogen ionization level. 
Similar Iv. iPauldrachI (|l987f) recognized that ionization from 
excited states shifts the ionization level of elements h igher 
in the wind of C Pup (04f), and iHillier et al] (|1983D and 
iHillieij (|l987t ) found that ionization from the n = 2 state of 
He was crucial in maintaining the He ionization balance in 
WR winds. 

During a CMFGEN or ASTAROTH simulation, the ra- 
diative recombination to and photo-ionization from all in- 
cluded levels can be monitored, so their importance can be 
easily assessed. Figure [5] shows the ratio of the total photo- 
ionization from the excited states of hydrogen to the total 
recombination to the same levels for both CMFGEN and 
ASTAROTH models. For the nebular approximation to be 
valid, the ratios in Fig. [5] need to be very small. However, 
the figure shows the ratios are significant, especially for 
the equatorial regions; therefore, we can conclude that the 
nebular approximation is not valid in the inner envelopes 
of these models. Ionization from the first few excited levels 
is extremely important, especially in the disk. 

To simulate the nebular approximation with 
ASTAROTH is not trivial. One needs to shut down 
the photo-ionization terms from excited states in the 
less dense regions, but needs to keep them in the hy- 
drostatic photosphere to recover LTE. Nevertheless, we 
experimented with spherical ASTAROTH models of the 
equatorial region of model A. To do this, we had to alter 
a few features of ASTAROTH. First, the inner boundary 
condition was changed from the diffusion approximation 
to a black-body function with Teff= 22,500 K. With 
this modification the inner hydrostatic envelope could 
be omitted and the radiative transition rates could be 
"safely" modified to simulate the nebular approximation. 
When we ran the above simulations (equatorial region of 
model A) with this version of the code, hydrogen became 
neutral down to ^1.1 R* with a negligible amount of HII. 
These experiments, therefore, fully supported the notion 
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Fig. 9. The neutral hydrogen fraction as a function of ra- 
dius in a spherical CMFGEN model created from the equa- 
torial region of model A (dashed line); and in the actual 
equatorial region of our model A (solid line). 



Despite the fact that hydrogen is ionized in most of 
our models, our simulations also prove that predominantly 
neutral material can exist around sgB[e] stars. In model 
D, for example, hydrogen has already recombined in the 
equatorial region. Furthermore, models B and C demon- 
strated that other species, like helium, can be neutral near 
the stellar surface even if hydrogen is ionized. It is a re- 
alistic expectation, therefore, that conditions suitable for 
molecule and dust formation may be found by more real- 
istic ASTAROTH models. The task for follow-up studies 
is to map the parameter space (much more extended than 
the one for spherical models) and find their location. Our 
calculations also highlighted the complex nature of sgB[e] 
envelopes and the limitations of simple semi-analytic cal- 
culations. For detailed mapping of the parameter space, 
self-consistent numerical simulations are necessary. 

Before we discuss our plans for the future, a few words 
about improving ASTAROTH and the sgB[e] models are 
warranted. One simplification in our sgB[e] models is the 
Sobolev approximation for the line transfer. At low lati- 
tudes the flow velocities were close to the ion thermal speed 
where the Sobolev approximation may no longer be accu- 
rate. We checked our results against CMFGEN models cre- 
ated from selected latitudes of the 2D models (e.g., pole 
or equator) that were using full line transfer. We saw no 
indication in these tests that departure from the Sobolev 
approximation would give a fundamentally different result. 
Note, that ASTAROTH and CMFGEN are very different 
codes and they use different solution techniques! 

Fig.[9]shows that the neutral hydrogen fraction is low in 
both the spherical non-Sobolev model and in the equator 
of the actual 2D calculation (model A) . Such a comparison 
is not entirely valid, but it can reveal large inconsistencies 
if present. The resuhs of CMFGEN and ASTAROTH, on 
the other hand, were in general agreement. The ionization 
level is approximately an order of magnitude lower in 2D 
since the radiation scatters off the dense disk preferentially 
toward the polar directions. This comparison also illustrates 
why modeling 2D problems with spherical codes, latitude 



Another simplification is that the effects of rotation 
were neglected in the simulations presented here. This 
is a significant omissio n if s gB[e] stars are fast rota- 
tors, as thought. iKrausI (|2006l ) found important changes 
when they introduced t hese effects in their earlier models 
()Kraus &: LamersI l2003t ) and we expect the same for our 
simulations. Rotation has four effects: First, rotation af- 
fects the dynamics of the circumstellar envelope, and the 
resultant line profiles. Second, rapid rotation causes gravity 
darkening su ch that the star i s hotter at the poles than at 
the equator (jvon ZeiDe]|[l92l . Third, very rapid rotation 
causes the star to become oblate. Fourth, rotation alters 
the escape probabilities and thus the occupation numbers 
of various levels. The inclusion of the fourth and the first 
two effects into ASTAROTH is easy; the third is difficult. 

We also anticipate significant changes when improving 
the atomic models. Iron is an especially important element 
to include because line-blanketing can seriously alter the 
radiation field, and hence influence the ionization of hy- 
drogen and helium. Unfortunately, this is a time consum- 
ing task because iron has tens-of-thousands of lines, and 
in order to allow for blanketing the line transfer must be 
done in the CMF. An alternative is to introduce the effects 
of line-blanke ting in an appro ximate manner (e.g., like in 
FASTWIND, iPuls et al.l[2005h . 

It is also very important to explore various dynamical 
concepts for sgB [e] envelopes and to include the non-radial 
velocity field. As mentioned in ijl] the bi-stability mech- 
anism is only one explanation for the observed bi-modal 
structure of the envelope. An important alternative is the 
Keplerian disk model which gained so me support by re- 
cent hydrodynamical simulations fe.g.. iMadura fc Owockil 
|2007[ ). The problem with models based on the classical line 
driven wind theory (jCastor et al.lll975[ ). like the bi-stability 
model, is that accelerating the large mass in the equato- 
rial regions requires multiple scattering which is difficult 
to maintain if the photons can easily escape in the po- 
lar directions. This is the very same effect that allowed 
lower hydrogen ionization in 2D than in the corresponding 
spherical models (see Fig.lHl). Furthermore, in hydrodynam- 
ical si mulation of sta.rs tha t are not close to the Eddington 
limit (jOwocki et al.lflQQSl ). the enhanced flux at the pole 
of a rapidly rotating star tends to create a bipolar outflow 
rather than an equatorial disk. 

In follow-up studies, we intend to improve ASTAROTH 
and will perform spectral analysis of specific sgB[e] stars. 
First, we plan to introduce non-radial velocities and the 
effect of rotation. We will also use a more realistic atomic 
model. The inclusion of iron blanketing is not trivial (pri- 
marily because of computational effort), so we will first 
concentrate on a realistic H, He, C, N, and O models. With 
these improvements we will undertake a spectroscopic anal- 
ysis of R 126, and will use the well-observed sarnple of 
LMC/SMC sgB[e] stars (Zic kgraf et al.lll985l Il986l 11991 
1996) to fully explore the allowed parameter ranges of the 
different sgB[e] models (Keplerian disk, bi-stability, etc) 
and look for observable differences. 

8. Summary 

Wc performed numerical calculations of hydrogen and he- 
lium ionization in sgBfel envelopes with ASTAROTH, our 
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of the sgB[c] parameter range. We also presented theoretical 
profiles for Ha and the C IV AA1550 doublet and examined 
their behavior for different viewing angles, and as a func- 
tion of the adopted stellar and envelope parameters. Our 
work yielded the following results: 

1. We found that it is much harder to form neutral hydro- 
gen in these sgB[e] envelopes than suggested by previ- 
ous studies. The reason for the discrepancy is the neb- 
ular approximation used by the earlier studies; our cal- 
culation shows that ionization from excited states of 
hydrogen is very important in the envelopes of sgB[e] 
stars and highlights the need for self-consistent numer- 
ical simulations to study these objects. 

2. Of the four models considered, only one formed a neu- 
tral hydrogen disk. Even in this case the disk did not 
extend to the stellar surface — it was truncated at ^2- 
3 R^,. A neutral hydrogen disk formed for the combi- 
nation of high mass-loss (> 10^''' M© yr"""^) and low 
effective temperature (Tg// < 18,000 K). 

3. An equatorial disk predominantly neutral in helium 
forms for all but the lowest mass-loss rates (< 
IQ-s Mq yr-i). 

4. The polar regions remain highly-ionized (e.g., CIV) in 
all models. 

5. The simulated spectra showed many familiar features 

of sgB[e] observations, like the narrow double-peak Ha 
emission and UV resonance line P Cygni profiles with- 
out emission. We found that Ha forms mainly in the 
disk, and the strength of the continuum scales with the 
projected area of the disk-wind interface. 

In the future we will include the effects of rotation (ini- 
tially gravity darkening, non-radial velocities; later surface 
distortions) and improve our atomic model by including 
N, O, and more levels of C. We will further explore the 
parameter range, allowed for sgB[e] stars, by using the 
well-observed LMC/SMC sample of these stars. An obvi- 
ous question to be addressed by these studies is - what is 
the structure of the disk (Keplerian, or dense outflow)? 
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